Cite this Article
Zhang Fu-Jie, Zhou Bao-Hua, Liu Xiao, Song Yu, Zuo Xu. Molecular dynamics simulation of atomic hydrogen diffusion in strained amorphous silica. Chinese Physics B, 2020, 29(2): 027101  
Permissions
Molecular dynamics simulation of atomic hydrogen diffusion in strained amorphous silica
Zhang Fu-Jie1, Zhou Bao-Hua1, Liu Xiao1, Song Yu2, 3, Zuo Xu1, 4, †
College of Electronic Information and Optical Engineering, Nankai University, Tianjin 300071, China
Microsystem and Terahertz Research Center, China Academy of Engineering Physics, Chengdu 610200, China
Institute of Electronic Engineering, China Academy of Engineering Physics, Mianyang 621999, China
Key Laboratory of Photoelectronic Thin Film Devices and Technology of Tianjin, Tianjin 300350, China

 

† Corresponding author. E-mail: xzuo@nankai.edu.cn

Project supported by the Science Challenge Project, China (Grant No. TZ2016003-1-105) and the CAEP Microsystem and THz Science and Technology Foundation, China (Grant No. CAEPMT201501).

Abstract

Understanding hydrogen diffusion in amorphous SiO2 (a-SiO2), especially under strain, is of prominent importance for improving the reliability of semiconducting devices, such as metal–oxide–semiconductor field effect transistors. In this work, the diffusion of hydrogen atom in a-SiO2 under strain is simulated by using molecular dynamics (MD) with the ReaxFF force field. A defect-free a-SiO2 atomic model, of which the local structure parameters accord well with the experimental results, is established. Strain is applied by using the uniaxial tensile method, and the values of maximum strain, ultimate strength, and Young’s modulus of the a-SiO2 model under different tensile rates are calculated. The diffusion of hydrogen atom is simulated by MD with the ReaxFF, and its pathway is identified to be a series of hops among local energy minima. Moreover, the calculated diffusivity and activation energy show their dependence on strain. The diffusivity is substantially enhanced by the tensile strain at a low temperature (below 500 K), but reduced at a high temperature (above 500 K). The activation energy decreases as strain increases. Our research shows that the tensile strain can have an influence on hydrogen transportation in a-SiO2, which may be utilized to improve the reliability of semiconducting devices.

PACS: 71.15.Pd;71.20.-b;66.30.-h;61.72.Hh
Keyword:molecular dynamics;tensile strain;amorphous SiO2;hydrogen diffusion
1. Introduction

The transportation of mobile impurities in dielectrics of semiconductor devices, such as metal–oxide–semiconductor field effect transistor (MOSFET), has been an intriguing topic, because it is closely related to the device reliability in harsh environments. Hydrogen is the most widely studied mobile species in amorphous silicon dioxide (a-SiO2), the most common material used as the gate dielectric in MOSFETs. Protons (positively charged hydrogen atoms) may drift toward Si/SiO2 interface in a manner of hopping between the fixed defects in a-SiO2 under a bias voltage. They will then de-passivate the hydrogen-saturated dangling bonds at the interface,[1,2] which generates the electrically-active dangling bond defects and results in device performance degrading.[3–5] Negative bias temperature instability (NBTI) is known as the most prevalent aging mechanism in MOSFETs, which increases the threshold voltage and limits the lifetime.[6–8] Reaction diffusion models for both atomic and molecular hydrogen have been proposed and verified experimentally to explain NBTI.[6,9] Strain is suspected to influence NBTI by distorting a-SiO2 lattice, and thus influencing hydrogen diffusion. However, the hydrogen atom diffusion in strained a-SiO2 is rarely investigated.

The experiments and theories mostly focused on hydrogen diffusion in a-SiO2 without strain. Using ab initio density-functional calculation, the energetics and dynamics of neutral hydrogen in α-crystobalite were explored.[10] It was discovered that a neutral hydrogen atom migrates by hopping between the local energy minima in the open voids. In addition, it was derived that the diffusivity is 8.1 × 10−3 cm2/s at 600 K and that the activation energy is 0.2 eV. The proton mobility was experimentally addressed by directly measuring the charge displacement,[11] and a short-time behavior involving an activation energy of 0.38 eV was discovered. Despite the limited range of studied temperature, the experiment strongly supported the calculation on the proton diffusion in a-SiO2,[12] thus concluding that the cross-ring inter-oxygen hopping assisted by network vibrations is the dominant diffusion mechanism of proton in a-SiO2 and the activation energy is 0.5 eV. Recently, by using MD simulations, the proton diffusion in a-SiO2 under strain was investigated in two temperature ranges.[6] It was shown that the activation energy is not influenced by the strain and varies around 0.07 ± 0.02 eV at low-temperature, and that on the contrary, the activation energy increases linearly with strain increasing at high-temperature. This result implied that the strain can be utilized to reduce proton diffusion in an a-SiO2 gate dielectric material. However, the mechanism behind the obvious discrepancy between the experimental values of hydrogen diffusivity was unclear, as the charge state of hydrogen that may affect the transportation was undetermined. In fact, it was shown by ab initio calculations that neutral hydrogen atom was thermodynamically unstable in a-SiO2.[13] However, the neutral hydrogen atom was experimentally detected in a-SiO2.[14] It was shown by ab-initio calculation that hydrogen atom in a-SiO2 can be neutral and non-bonded interstitial atom, positively charged and bonded to O atom, or negatively charged and bonded to Si atom.[15,16]

In this work, hydrogen diffusion in a-SiO2 under strain is simulated by using classical molecular dynamics (MD). A diffusion mechanism parallel to the proton hopping is proposed for the neutral hydrogen atoms in a-SiO2. The diffusivity and activation energy of hydrogen atom in strained a-SiO2 are derived from the transportation paths simulated by MD. It is shown that applying 3% tensile strain can remarkably reduce the diffusivity at high temperature (above 500 K), although it lowers the activation energy to 1/3 of the value in the case of no strain. This effect may be utilized to control hydrogen diffusion in semiconducting devices.

2. Model and methods
2.1. a-SiO2 model

The a-SiO2 model was first constructed by simulating the melting and subsequently quenching of crystalline silica (c-SiO2) through using the Large-scale atomic/molecular massively parallel simulator[17] (LAMMPS) with the ReaxFF force field.[18] The unit cell of cristobalite silica[19,20] with 1152 atoms and a size of 20.25 Å × 30.50 Å × 28.15 Å was first heated at 8000 K for 200 ps with an isothermal and isochoric (NVT) ensemble. The heating temperature and duration were sufficient to melt the c-SiO2 unit cell and completely remove the memory of the initial crystalline structure. The unit cell was subsequently quenched from 8000 K to room temperature (300 K) by using an NVT ensemble at 5 K/ps.[21–25] The cooled silica was further relaxed at 300 K and 1 atm for 200 ps with an isothermal and isobaric (NPT) ensemble to approach to the equilibrium structure of a-SiO2. Nose–Hoover thermostat and Berendsen barostat were used to control the temperature and pressure in the simulations. The MD time step was 0.5 fs in all simulations. The a-SiO2 structure generated by classical MD was then optimized at 0 K by using ab-initio calculation based on the density functional theory (DFT). The structure relaxation was implemented by the Vienna ab initio simulation package (VASP).[26] The cut-off energy was set to be 400 eV, and breaking conditions of electronic and ionic loops were 10−3 eV and 10−2 eV/Å, respectively. The Brillouin-zone (BZ) integration was performed only at the Γ-point, since the a-SiO2 unit cell was sufficiently large. The initial c-SiO2 and the final a-SiO2 structures are shown in Fig. 1.

Fig. 1. Atomistic model of (a) c-SiO2 and (b) a-SiO2 containing 1152 atoms.
2.2. Uniaxial strain

The a-SiO2 model was loaded under uniaxial tension along the x-axis. The strain ε was introduced by stretching the a-SiO2 model, and defined as where a and a0 are the loaded and unloaded lattice parameter, respectively. The uniaxial tensions were simulated at 300 K with an NVT ensemble and parallelepiped periodic boundary conditions were adopted along all three directions.[18,25] During the MD simulation, the lattice parameters along the two directions perpendicular to the strain were fixed at their initial values, respectively. Four different stress–strain curves were obtained at four strain rates from 1 × 1011 s−1 to 1 × 1014 s−1, and the elastic properties (including maximum strain, ultimate strength, and Young’s modulus) were derived from these curves. All simulations were carried out with the identical setups to enable a systematic comparison.

2.3. Diffusivity and activation energy

The diffusivity was calculated from the mean square displacement (MSD) derived from the trajectories of seven hydrogen atoms randomly placed in the a-SO2 sample.[6,27–29] The simulations were carried out at four different temperatures (500 K, 800 K, 1000 K, 1500 K) using an NVT ensemble with a duration of 300 ps in time steps of 0.5 fs. The MSD was computed from the trajectories of the hydrogen atoms traced in the simulations by the equation where the angle brackets refer to the ensemble average implemented by averaging over the truncated trajectories with a time span t0. The diffusion coefficient D at temperature T was then derived from the MSD following the Einstein relationship,

The calculated D(T) was then fitted to the Arrhenius form[30]
where D0 and Ea are the pre-exponential factor and the activation energy, respectively.

3. Results and discussion
3.1. a-SiO2 model

There is no coordination defect in the a-SiO2 model, i.e., each Si atom is bonded to four O atoms and each O atom is bonded to two Si atoms. The density, average bond length and angle, and the peaks of the radial distribution function (RDF) of the a-SiO2 model are listed in Table 1, which are consistent well with the measurements and the MD simulations. It is worth noting that the density of the a-SiO2 model is slightly lower than the experimental value, which can be caused by the cooling rate.[31–33] The first O–O peak is 2.89% smaller than the experimental value, and the first Si–Si and Si–O peak position are 2.37% and 2.61% larger than the experimental values, respectively. These discrepancies can be attributed to the statistical fluctuation due to the relatively small size of the model.

Table 1.

Structural properties of a-SiO2 model.

.
3.2. Strain effect

The stress–strain curves of the a-SiO2 model under the uniaxial tension are shown in Fig. 2. Four different strain rates are adopted in the calculation to compare with previously reported results.[18,25,38–41] The maximum strain, where the uniaxial tensile simulation stops, is 0.8, far exceeding the strain at the maximum stress. By increasing the tensile rate from 1 × 1011 s−1 to 1 × 1014 s−1, the strain value at the maximum stress is enhanced from 21.7% to 48.9%, and the ultimate strength increases from 25.28 GPa to 39.02 GPa. This indicates that the a-SiO2 model exhibits a higher strength and strain to failure as the strain rate increases.

Fig. 2. Effects of strain rate on stress–strain response of a-SiO2 model.

The stress–strain curves show a linear elastic behavior if the strain is less than 1% (Fig. 2), and the Young’s modulus is derived (Table 2). In addition, the strain and ultimate strength at the maximum stress for different strain rates are listed in Table 2. The maximum strain is 21.7% at a strain rate of 1 × 1011 s−1, close to the experimental value (23%).[41,42] In addition, the Young’s modulus matches the experimental value of silicate glass (71.9 GPa).[41,42] Therefore, we select the strain rate of 1 × 1011 s−1 for further study.

Table 2.

Mechanical properties of a-SiO2 model at different strain rates.

.

The influences of the tensile strain on the Si–O bond length, and the O–Si–O and Si–O–Si bond angle under the strain values of 0%, 5%, 10%, 15%, and 20% are displayed in Fig. 3, where the Si–O bond length and the Si–O–Si bond angle show a sensitive dependence on uniaxial tension. With the strain increasing the Si–O bond length increases from 1.61 Å to 1.65 Å, and the Si–O–Si angle increases from 157° to 164°, as shown in Figs. 3(a) and 3(b). The O–Si–O angle is stiffer than the Si–O–Si angle, and its value is kept around 109° in Fig. 3(c). These results are in good agreement with the calculations.[18,42,43]

Fig. 3. Influence of the tensile strain on local structures of a-SiO2 model: variations of probability with (a) Si–O bond length, (b) Si–O–Si bond angle, and (c) O–Si–O bond angle.
3.3. Hydrogen atom diffusion in unstrained a-SiO2

The diffusion trajectory of a single hydrogen atom at 1000 K is shown in Fig. 4(a). The hydrogen atom initially vibrates at an equilibrium location for more than 200 ps, then hops to another equilibrium location, and keeps vibrating around that location until the next hop.

Fig. 4. (a) Hydrogen atom diffusion over a time span of 300 ps, and (b) hydrogen atom at local energy minimum and its means distances to the neighboring oxygen atoms, with pink, cyan, and purple spheres denoting Si, O, and H atoms, respectively.

It was shown by the previous calculations that neutral hydrogen atom is thermodynamically unstable in a-SiO2 and that either positively or negatively charged hydrogen atom is stable, depending on the Fermi level.[13] In the DFT simulation of proton diffusion, the proton was observed to vibrate around an equilibrium position and the trajectory of the proton vibration was approximately in the Si–O–Si plane with an average H–O distance of about 1.1 Å.[12] The proton hopping was assisted by the vibration of the Si–O framework. The hopping to the nearest-neighbor oxygen was infrequent, compared with the crossing-ring hopping, in which process the shortest O–O distance was about 2.5 Å. The hopping mechanism was confirmed by using classical MD simulations,[6] and it was concluded that the hydrogen atom diffused in a-SiO2 in the form of the proton. However, it was observed that the average H–O length was 2.1 Å and the crossing-ring O–O distance involved in the hopping was about 4.2 Å in the simulations.

Our simulations, however, show that the proton hopping mechanism can be paralleled by another mechanism, where the hydrogen atom is in neutral charge state. In fact, the hydrogen atom will be neutral, if located at the local energy minimum in an oxygen cage, where it keeps vibrating for a long time (several hundred picoseconds) but is not bonded to any oxygen atom (Fig. 4(b)). The same behavior was previously discovered in the study of the diffusion of neutral hydrogen in α-crystobalite.[10] This local energy minimum of the hydrogen atom is similar to that near the tetrahedral interstitial site in crystalline Si due to the interaction between the hydrogen atom and the Si lattice. The spin charge density is illustrated in Fig. 5, where it is distributed around the hydrogen nucleus, implying that there is an electron on the hydrogen atom during diffusion. The charge state of the hydrogen atom is further identified by using Bader charge analysis. The Bader charge of the hydrogen atom is |0.99|. In addition, it is observed that the O–H distance is much longer than 1.1 Å, the average O–H distance for a proton bonding an oxygen bridge. The distances between the hydrogen atom and the surrounding oxygen atoms are marked in Fig. 6, where the O–H distance is around 2.0 Å in the vibration, far exceeding the H–O bond length of 1.1 Å. Therefore, we predict that there should be an alternative diffusion mechanism paralleling to the cross-ring proton hopping mechanism, where the hydrogen atom hops from a local energy minimum in an oxygen cage like that in Fig. 4(b) to an adjacent one as indicated by the O–H distances in Fig. 6. It can be observed in Fig. 6 that the hydrogen atom vibrates at a local energy minimum in a time from 0 to 225 ps and then hops to an adjacent energy minimum.

Fig. 5. Spin charge density of hydrogen atom diffusing in the a-SiO2 network, where the isosurface value is set to be 0.01 e/Å3, with pink, cyan, and purple spheres denoting Si, O, and H atoms, respectively.
Fig. 6. Evolution of H–O distance in 300 ps simulation of hydrogen atom diffusion, where the subscripted oxygen atoms are illustrated in Fig. 4(b).

Diffusivity as a kinetic quantity can be written in Arrhenius form, which is determined by two important parameters: activation energy and pre-exponential factor (D0). The activation energy measured in experiment spreaded from 0.05 eV to 0.9 eV.[6,10–12,44,45] Specifically, for neutral atomic hydrogen, in dense oxides the activation energy was measured to be between 0.1 eV and 0.2 eV, whereas in open silica channels the barrier was found to be lower, about 0.05 eV.[10,14,46,47] However, the pre-exponential factor is experimentally measured rarely.[10] Specifically, in wet fused silica, the pre-exponential of hydrogen atom was estimated at 1 × 10−4 cm2/s in theory.[10,14,33] The logarithm of the calculated diffusivity (ln D) is plotted and linearly fitted against inverse temperature (1/kT) in Fig. 7, from which the activation energy of 0.57 eV and the pre-exponential factor of 5.5 × 10−4 cm2/s are derived. The range of activation energy and pre-exponential factor are in agreement with our calculated values.

Fig. 7. Arrhenius plot of hydrogen atom diffusivity in a-SiO2.
3.4. Hydrogen atom diffusion in strained a-SiO2

The influences of the strain on the calculated diffusivity at different temperatures is shown in Fig. 8, where the diffusivity values are calculated respectively at the tensile strains of 3%, 5%, 7%, and 10%. At a relatively low temperature (500 K), the diffusivity shows a general trend of increase with strain increasing, except the decrease in a range from 3% to 5% strain. More than that, the diffusivity increases by 9.7 times at 3% strain, which implies that a moderate tensile strain may significantly enhance the hydrogen atom diffusion at a relatively low temperature. On the contrary, the simulated diffusivity shows a general trend of decrease with strain increasing at the simulation temperatures higher than 500 K. In addition, the diffusivity drastically decreases as the strain increases from 0 to 3%. This diffusivity is about 4.6 times less than that at 800 K, and about one order of magnitude less than that at 1000 K and 1500 K. With strain further increasing beyond 3%, the diffusivity fluctuates at 800 K and becomes almost flat at 1000 K and 1500 K.

Fig. 8. Dependence of the calculated hydrogen diffusivity on strain at different temperatures.

The calculated activation energy of hydrogen atom diffusion is plotted against strain in Fig. 9(a), showing a decreasing tendency with strain increasing. It drastically decreases from 0.57 eV to 0.22 eV, about 1/3, as the strain increases from 0 to 3%. This result is consistent with that from the model of neutral hydrogen atom diffusion. In that model, the hydrogen atom vibrates around a local energy minimum bound in an oxygen atom cage. The diffusion is microscopically accomplished by a series of hops from a cage to an adjacent one. Each hop needs to pass through an Si–O ring perpendicularly and overcomes an energy barrier associated with the activation energy statistically. When a tensile strain is applied, the Si–O bond length and the Si–O–Si bond angle increase on average, statistically implying that the Si–O rings enlarge and thus the energy barrier associated with the hydrogen atom hopping decreases.

Fig. 9. (a) Calculated activation energy of hydrogen atom diffusion as a function of applied strain, and (b) reaction curve of hydrogen atom diffusion between two adjacent local energy minima at different values of strain.

The strain effect on the activation energy is further investigated by calculating the microscopic diffusion barriers at different strains in Fig. 9(b). The reaction curves are calculated for a diffusion event at four different strains by using the climbing image nudged elastic band (CI-NEB) method with a spring constant of 5.0 eV/Å between adjacent images.[48] For each curve, the initial and final states of the curves are two local energy minima located in two adjacent oxygen cages, respectively. The transition state is corresponding to the structure where the hydrogen atom penetrates through the Si–O ring separating the two cages, or mathematically the saddle point between the two local energy minima. The forward (or reverse) energy barrier is calculated to be the difference between the initial (or final) state and the transition state.[49] As the strain increases from 0% to 5%, the forward energy barrier monotonically decreases from 0.43 eV to 0.36 eV. The energy barriers derived from the reaction curves are in the same order of magnitude as the activation energy values derived from the diffusion trajectories. This is qualitatively consistent with the decreasing trend of the activation energy with strain increasing.

4. Conclusions

We perform atomic-scale simulations to investigate the diffusion of hydrogen atom in a-SiO2 and the influence of strain on diffusivity and activation energy. The a-SiO2 model is prepared and its structural parameters (mass density, average Si–O bond length, Si–O–Si and O–Si–O bond angles, and the first-peaks of RDFs) are compared with the measurements and calculations. Uniaxial strain is then applied to the model at different tensile rates. The stress–strain response stiffens and achieves higher strength and strain to failure as the tensile rate increases. In addition, we also find Young’s modulus enhancement with the increase of tensile rate. The diffusions of hydrogen atom in a-SiO2 are simulated at different values of strain. The simulations reveal another diffusion picture paralleling to the well-known cross-ring proton hopping, that is, the hydrogen atom mostly remains in neutral charge state and vibrates at a local energy minimum in an oxygen cage and hops to another local energy minimum in the adjacent cage assisted by the vibration of the a-SiO2 framework. By fitting the simulated diffusivity to the Arrhenius plot, the diffusion activation energy of the hydrogen atom is derived to be 0.57 eV in the case of no strain and decreases as the strain increases. Specially, it drops to about 1/3 of the intimal value as the strain increases from 0 to 3%. In addition, the diffusivity shows opposite trend with increasing strain, depending on the temperature. It generally increases as the strain increases at 500 K but decreases at 1000 K and 1500 K. More than that, applying 3% strain can drastically lower the diffusivity down by about 13.3 times at 1500 K. This implies that the strain may be utilized to efficiently inhibit hydrogen diffusion at a relatively high temperature. This research may help to understand hydrogen atom diffusion in a-SiO2 on an atomic scale and to optimize the fundamental processing of semiconducting devices based on Si/a-SiO2 interface, for which the reliability issues are directly related to hydrogen and its diffusion in a-SiO2 are critical.

Reference
[1] Sheikholeslam S A Manzano H Grecu C Ivanov A 2018 Superlattices Microstruct. 120 561
[2] Pantelides S T Rashkeev S N Buczko R Fleetwood D M Schrimpf R D 2000 IEEE Trans. Nucl. Sci. 47 2262
[3] Witczak S C Suehle J S Gaitan M 1992 Solid State Electron. 35 345
[4] Vanheusden K Devine R A B 2000 Appl. Phys. Lett. 76 3109
[5] Afanas’ev V V Adriaenssens G J Stesmans A 2001 Microelectron. Eng. 59 85
[6] Sheikholeslam S A Manzano H Grecu C Ivanov A 2016 J. Mater. Chem. 4 8104
[7] Alam M A 2003 IEEE Int. Electron Devices Meet. December 8–10, 2003 Washington DC USA 14.4.1 https://doi.org/10.1109/IEDM.2003.1269295
[8] Krishnan A T Chancellor C Chakravarthi S Nicollian P E Reddy V Varghese A Khamankar R B Krishnan S 2005 IEEE Int. Electron Devices Meet. Tech. Dig. December 5, 2005 Washington, DC, USA 688 https://doi.org/10.1109/IEDM.2005.1609445
[9] KflogluKufluoglu H Alam M A 2007 IEEE Trans. Electron. Dev. 54 1101
[10] Tuttle B 2000 Phys. Rev. 61 4417
[11] Devine R A B Herrera G V 2001 Phys. Rev. 63 233406
[12] Godet J Pasquarello A 2006 Phys. Rev. Lett. 97 155901
[13] Godet J Pasquarello A 2005 Microelectron. Eng. 80 288
[14] Griscom D L 1985 J. Appl. Phys. 58 2524
[15] Yue Y Wang J Zhang Y Song Y Zuo X 2018 Physica 533 5
[16] Xiong K Robertson J Clark S J 2007 J. Appl. Phys. 102 083710
[17] Plimpton S 1995 J. Comput. Phys. 117 1
[18] Duan F L Zhang C Liu Q S 2015 J. Nano Res. 30 59
[19] Petousis I Mrdjenovich D Ballouz E Liu M Winston D Chen W Graf T Schladt T D Persson K A Prinz F B 2017 Sci. Data 4 160134
[20] Mathew K Zheng C Winston D Chen C Dozier A Rehr J J Ong S P Persson K A 2018 Sci. Data 5 180151
[21] Ebrahem F Bamer F Markert B 2018 Comput. Mater. Sci. 149 162
[22] Smedskjaer M M Bauchy M Mauro J C Rzoska S J Bockowski M 2015 J. Chem. Phys. 143 164505
[23] Wang B Yu Y Wang M Mauro J C Bauchy M 2016 Phys. Rev. 93 064202
[24] Yu Y Krishnan N M A Smedskjaer M M Sant G Bauchy M 2018 J. Chem. Phys. 148 074503
[25] Chowdhury S C Haque B Z Gillespie J W 2016 J. Mater. Sci. 51 10139
[26] Kresse G 2007 http://cms.mpi.univie.ac.at/vasp/
[27] Wood S M Eames C Kendrick E Islam M S 2015 J. Phys. Chem. 119 15935
[28] van Duin A C Merinov B V Han S S Dorso C O Goddard W A 2008 J. Phys. Chem. 112 11414 3rd
[29] Sheikholeslam S A Luo W Grecu C Xia G Ivanov A 2016 J. Non Cryst. Solids 440 7
[30] Bauer T Lunkenheimer P Loidl A 2013 Phys. Rev. Lett. 111 225702
[31] Sundararaman S Ching W Y Huang L 2016 J. Non Cryst. Solids 445�?46 102
[32] Lane J M 2015 Phys. Rev. 92 012320
[33] Vollmayr K Kob W Binder K 1996 Phys. Rev. 54 15808
[34] Mozzi R L Warren B E 1969 J. Appl. Crystallogr. 2 164
[35] Fogarty J C Aktulga H M Grama A Y van Duin A C Pandit S A 2010 J. Chem. Phys. 132 174704
[36] El-Sayed A M Watkins M B Afanas’ev V V Shluger A L 2014 Phys. Rev. 89 125201
[37] Da Silva J R G Pinatti D G Anderson C E Rudee M L 1975 Philos. Mag. 31 713
[38] Zhang Y Huang L Shi Y 2019 Nano Lett. 19 5222
[39] Yuan F Huang L 2012 J. Non Cryst. Solids 358 3481
[40] Muralidharan K Oh K D Deymier P A Runge K Simmons J H 2007 J. Mater. Sci. 42 4159
[41] Gupta P K Kurkjian C R 2005 J. Non Cryst. Solids 351 2324
[42] Muralidharan K Simmons J H Deymier P A Runge K 2005 J. Non Cryst. Solids 351 1532
[43] Pedone A Malavasi G Menziani M C Segre U Cormack A N 2008 Chem. Mater. 20 4356
[44] Beyer W Wagner H 1982 J. Appl. Phys. 53 8745
[45] Fink D Krauser J Nagengast D Murphy T A Erxmeier J Palmetshofer L Bräunig D Weidinger A 1995 Appl. Phys. 61 381
[46] Verdi L Miotello A 1993 Phys. Rev. 47 14187
[47] Cartier E Buchanan D A Stathis J H DiMaria D J 1995 J. Non Cryst. Solids 187 244
[48] Henkelman G Uberuaga B P Jónsson H 2000 J. Chem. Phys. 113 9901
[49] Yue Y Song Y Zuo X 2018 Chin. Phys. 27 037102